home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / pascal / p_mat.exe / PMAT.PAS < prev    next >
Pascal/Delphi Source File  |  1993-01-23  |  24KB  |  895 lines
  1. {
  2. P-Mat - A Turbo Pascal program for Recursive Matrix Algebra
  3.  
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/30/93
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9.  
  10.  
  11. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  12. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  13. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  14. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  15. FROM USE OF THIS PROGRAM.
  16. }
  17.  
  18. Unit pmat;
  19.  
  20. Interface
  21.  
  22. Type
  23.      fp = ^double;
  24.      ap = Array[1..1] Of double;
  25.      apptr = ^ap;
  26.      app = Array[1..1] Of apptr;
  27.  
  28.      vmatrixptr = ^vmatrix;
  29.      vmatrix = Object
  30.         r, c : integer;
  31.         Function m( i, j: integer ): double;
  32.         Function mm( i, j: integer ) : fp;
  33.         constructor MakeMatrix( vr, vc: integer );
  34.         Destructor KillVmatrix;
  35.         Procedure Garbage;
  36.         Procedure Show( strng: String );
  37.         Procedure infomatrix( strng: String );
  38.  
  39.         private
  40.            v,vcheck: ^app;
  41.            nelems : longint;
  42.            onstack : boolean;
  43.            Procedure purgevectors;
  44.            Procedure allocvectors( rr, cc: integer );
  45.      End;
  46.  
  47.      mstackptr = ^mstack;
  48.      mstack = Object
  49.         constructor InitMatrixStack;
  50.         Destructor KillMatrixStack;
  51.         Procedure inclevel;
  52.         Procedure declevel;
  53.         Function ReturnMat: vmatrixptr;
  54.         Function DecReturn: vmatrixptr;
  55.         Procedure push( Var a: vmatrixptr );
  56.         Procedure nrerror( strng: String );
  57.         Procedure DumpStack;
  58.  
  59.         private
  60.           nummats, level : integer;
  61.           next : mstackptr;
  62.           mv : vmatrixptr;
  63.           Procedure cleanstack( a: vmatrixptr );
  64.           Function pop: vmatrixptr;
  65.      End;
  66.  
  67. Var
  68.     dispatch : mstackptr;
  69.  
  70. Function matequals( Var lop: vmatrixptr; rop: vmatrixptr ) : vmatrixptr;
  71. Function ident( n: integer ): vmatrixptr;
  72.  
  73. Function Mult( A, B: vmatrixptr ):vmatrixptr;
  74. Function emult( a, b: vmatrixptr ):vmatrixptr;
  75. Function scmult( a: double; b: vmatrixptr ): vmatrixptr;
  76. Function multsc( b: vmatrixptr; a: double ): vmatrixptr;
  77.  
  78. Function divis( a, b: vmatrixptr ):vmatrixptr;
  79. Function scdivis( a: double; b: vmatrixptr ): vmatrixptr;
  80. Function divissc( b: vmatrixptr; a: double ): vmatrixptr;
  81.  
  82. Function add( a, b: vmatrixptr ): vmatrixptr;
  83. Function scadd( a: double; b: vmatrixptr ): vmatrixptr;
  84. Function addsc( b: vmatrixptr; a: double ): vmatrixptr;
  85.  
  86. Function neg( A: vmatrixptr ): vmatrixptr;
  87. Function sub( A, B: vmatrixptr ):vmatrixptr;
  88. Function scsub( A: double; B: vmatrixptr ):vmatrixptr;
  89. Function subsc( B: vmatrixptr; A: double ):vmatrixptr;
  90.  
  91. Function tran( a: vmatrixptr ): vmatrixptr;
  92. Function sweep( A: vmatrixptr; k1, k2: integer ) : vmatrixptr;
  93. Function inv( Amat: vmatrixptr ): vmatrixptr;
  94.  
  95. Function fill( rr, cc: integer; d: double ):vmatrixptr;
  96. Function submat( a: vmatrixptr; lr, r, lc, c: integer ): vmatrixptr;
  97. Function cv( a, b: vmatrixptr ): vmatrixptr;
  98. Function ch( a, b: vmatrixptr ): vmatrixptr;
  99. Function vecdiag( a: vmatrixptr ): vmatrixptr;
  100.  
  101. Function reada( fid: String ): vmatrixptr;
  102. Procedure writea( fid: String; a: vmatrixptr; titlestr: String );
  103.  
  104. Procedure setwid( wid: integer );
  105. Procedure setdec( decimals: integer );
  106. Function getwid: integer;
  107. Function getdec: integer;
  108.  
  109. Implementation
  110.  
  111. {///////////////// vmatrix functions ///////////////}
  112. {$R-}
  113. Procedure vmatrix.allocvectors;
  114. Var i: integer;
  115. Begin
  116.     if c > 8192 then 
  117.        dispatch^.nrerror('A matrix cannot have more than 8192 columns');
  118.     getmem( v, r * sizeof( apptr ) );
  119.     For i := 1 To r Do
  120.         getmem( v^[i], c * sizeof( double ) );
  121. End;
  122.  
  123. Procedure vmatrix.purgevectors;
  124. Var i: integer;
  125. Begin
  126.     If v = Nil Then exit;
  127.     For i := r Downto 1 Do
  128.         freemem( v^[i], c * sizeof( double ) );
  129.     freemem( v, r * sizeof( apptr ) );
  130. End;
  131. {$R+}
  132.  
  133. constructor vmatrix.MakeMatrix;
  134. Begin
  135.     r := vr;
  136.     c := vc;
  137.     onstack := false;
  138.     if c > 8192 then 
  139.        dispatch^.nrerror('A matrix cannot have more than 8192 columns');
  140.     nelems := longint( longint( r ) * longint( c ) ) * longint( sizeof( double ) );
  141.     If nelems < maxAvail - 256 Then allocvectors( r, c )
  142.     Else dispatch^.nrerror( 'cannot allocate matrix' );
  143.     vcheck := v;
  144. End; { MakeMatrix }
  145.  
  146. Destructor vmatrix.KillVmatrix;
  147. Begin
  148.     purgevectors;
  149.     vcheck := Nil;
  150. End;
  151.  
  152. Procedure vmatrix.Garbage;
  153. Var errorint : boolean;
  154. Begin
  155.     errorint := false;
  156.     If vcheck <> v Then errorint := true;
  157.     If v = Nil Then errorint := true;
  158.     If r < 1 Then errorint := true;
  159.     If c < 1 Then errorint := true;
  160.     If errorint Then dispatch^.Nrerror( 'garbage' );
  161. End;
  162.  
  163. {$R-}
  164.  
  165. Function vmatrix.m;
  166. Begin
  167.     { access a member on the right side of assignment }
  168.     If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
  169.         dispatch^.nrerror( 'm: index out of range' );
  170.     m := v^[i]^[j];
  171. End;
  172.  
  173. Function vmatrix.mm;
  174. Begin
  175.     { store a matrix element on the left side of assignment }
  176.     If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
  177.         dispatch^.nrerror( 'mm: index out of range' );
  178.     mm := @v^[i]^[j];
  179. End;
  180.  
  181. {$R+}
  182.  
  183. Procedure CopySD( Var source, dest: vmatrixptr );
  184. Var
  185.   i,j : integer;
  186. Begin
  187.     source^.Garbage;
  188.     dest^.Garbage;
  189.     
  190.     If source = dest Then exit;
  191.     If source^.onstack Then
  192.         With dest^ Do Begin 
  193.             purgevectors;
  194.             nelems := source^.nelems;
  195.             v := dispatch^.Pop^.v;
  196.             r := source^.r;
  197.             c := source^.c;
  198.             vcheck := v;
  199.             source^.v := Nil;
  200.             dispose( source, killvmatrix );
  201.             exit;
  202.         End;
  203.     
  204.     If source^.v = dest^.v Then
  205.         With dest^ Do Begin 
  206.             r := source^.r;
  207.             c := source^.c;
  208.             allocvectors( r, c );
  209.             vcheck := v;
  210.             nelems := source^.nelems;
  211.         End;
  212.     
  213.     If ( dest^.r <> source^.r ) Or( dest^.c <> source^.c ) Then
  214.         With dest^ Do Begin 
  215.             purgevectors;
  216.             r := source^.r;
  217.             c := source^.c;
  218.             allocvectors( r, c );
  219.             vcheck := v;
  220.             nelems := source^.nelems;
  221.         End;
  222.     
  223.     For i := 1 To dest^.r Do
  224.         For j := 1 To dest^.c Do
  225.             dest^.mm( i, j )^ := source^.m( i, j );
  226. End;
  227.  
  228. Function matequals;
  229.      Begin
  230.          rop^.Garbage;
  231.          lop^.Garbage;
  232.          copySD( rop, lop );
  233.          dispatch^.CleanStack( lop );
  234.          matequals := lop;
  235.      End;
  236.  
  237.  
  238. {////////////////////// stack functions /////////}
  239.  
  240.  
  241. constructor mstack.InitMatrixStack;
  242. Begin
  243.     nummats := 0;
  244.     level := 1;
  245.     mv := Nil;
  246.     getmem( next, sizeof( mstack ) );
  247.     With next^ Do Begin 
  248.         level := 0;
  249.         next := Nil;
  250.         mv := Nil;
  251.         nummats := 0;
  252.     End;
  253. End;
  254.  
  255. Destructor mstack.KillMatrixStack;
  256. Begin
  257.     freemem( next, sizeof( mstack ) );
  258. End;
  259.  
  260. Procedure mstack.nrerror;
  261. Begin
  262.     writeln( 'Fatal Error: ', strng );
  263.     writeln( 'Exiting to system' );
  264.     halt;
  265. End;
  266.  
  267. Procedure mstack.inclevel;
  268. Begin
  269.     level := level + 1;
  270. End;
  271.  
  272. Procedure mstack.declevel;
  273. Begin
  274.     level := level - 1;
  275.     If level <= 0 Then nrerror( 'declevel: decremented too often' );
  276. End;
  277.  
  278. Function mstack.ReturnMat;
  279. Begin
  280.     ReturnMat := next^.mv;
  281. End;
  282.  
  283. Function mstack.DecReturn;
  284. Begin
  285.     declevel;
  286.     DecReturn := next^.mv;
  287. End;
  288.  
  289. Procedure mstack.push;
  290.    Var
  291.       t : mstackptr;
  292. Begin
  293.     a^.Garbage;
  294.     getmem( t, sizeof( mstack ) );
  295.     t^.level := dispatch^.level;
  296.     t^.next := dispatch^.next;
  297.     t^.mv := a;
  298.     a^.onstack := true;
  299.     dispatch^.next := t;
  300.     dispatch^.nummats := dispatch^.nummats + 1;
  301.     t^.nummats := dispatch^.nummats;
  302. End;
  303.  
  304. Function mstack.pop;
  305. Var t : mstackptr;
  306.     vv : vmatrixptr;
  307. Begin
  308.     If dispatch^.level < 0 Then
  309.         nrerror( 'pop: missmatched inclevel-declevel statments, level < 0' );
  310.     If ( dispatch^.next^.next = Nil ) Or( Nil = dispatch^.next^.mv ) Then
  311.         nrerror( 'pop: popping off the tail' );
  312.     
  313.     t := dispatch^.next;
  314.     dispatch^.next := di